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ABSTRACT 

The long-term evolution of stellar orbits bound to a massive centre is studied in 
order to understand the cores of star clusters in central regions of galaxies. Stellar 
trajectories undergo tiny perturbation, the origin of which is twofold: (i) gravitational 
field of a thin gaseous disc surrounding the galactic centre, and (ii) cumulative drag 
due to successive interactions of the stars with material of the disc. Both effects are 
closely related because they depend on the total mass of the disc, assumed to be a 
small fraction of the central mass. It is shown that, in contrast to previous works, 
most of the retrograde (with respect to the disc) orbits are captured by the central 
object, presumably a massive black hole. Initially prograde orbits are also affected, so 
that statistical properties of the central star cluster in quasi-equilibrium may differ 
significantly from those deduced in previous analyses. 

Key vifords: accretion, accretion discs - galaxies: nuclei - celestial mechanics, stellar 
dynamics 



1 INTRODUCTION 

This paper extends previous studies on interaction between 
stars and an accretion disc near a massive galactic nucleus. 
Relevant references are, in particular, Syer, Clarke & Rees 
(1991, these authors estimate time-scales for the evolution 
of stellar orbital parameters in the Newtonian regime), and 
Vokrouhlicky & Karas (1993, relativistic generalization deal- 
ing with individual trajectories). Pineault & Landry (1994) 
and Rauch (1995) studied statistical properties of stellar or- 
bits in a dense cluster near a galactic core with an accretion 
disc. 

Observational evidence and theoretical considerations 
suggest that many galaxies harbour very massive compact 
cores (Mc ~ 10^-10^ M0), presumably black holes. In partic- 
ular, high energy output, variability, spectral properties, and 
production of jets in active galactic nuclei (AGN) can be un- 
derstood in terms of the model with a supermassive central 
object surrounded by an accretion disc (e.g., Courvoisier & 
Mayor 1990; Urry & Padovani 1995). However, linear resolu- 
tion of present observational techniques corresponds at best 
to several hundreds of gravitational radii of the hypothetical 
black hole. The innermost regions of these galaxies cannot 
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thus be resolved and conclusions about their structure must 
be inferred from integral characteristics (integrated over the 
angular and temporal resolution of the device used for ob- 
servation). Distribution of stars and gaseous material close 
to the galactic centre is one of the important tools in this 
respect because velocity dispersion and the corresponding 
luminosity profile of the nucleus reflect the presence and 
properties of the central massive object and the disc (Perry 
& Williams 1993; cf. Marconi et al. 1997 for recent observa- 
tional results). 

We will study the situation in which the central object 
is surrounded by an accretion disc and a dense star clus- 
ter. It is the aim of the present contribution to examine the 
role of periodic interactions of the stars with the disc mate- 
rial, simultaneously considering the gravitation of the disc. 
Mutual gravitational interaction of stars forming a dense 
cluster has been studied since the early works of Ambart- 
sumian (1938) and Spitzer (1940) while the importance of 
star-disc collisions for the structure of galactic nuclei has 
been recognized in the early 1980s (Goldreich & Tremaine 
1980; Ostriker 1983; Hagio 1987). Huang & Carlberg (1997) 
studied a related problem in the dynamics of galaxies. 

The gravitation of accretion discs was neglected in pre- 
vious works because its mass, Md, is presumed very small 
compared to the mass of the central object (fi = M^/Mc <^ 
1; ^ is a free parameter in our study). We also assume ^ ^ 1 
so that the gravity of the disc acts as a perturbation on the 
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stellar motion around the central mass. We will show, inde- 
pendent of the precise value of ^, that the effect of the disc 
gravity on the long-term evolution of stellar orbits must be 
taken into account together with star-disc interactions. In 
particular, we will show that circularization of many of the 
orbits, evolution of their inclination, and stellar capture rate 
are visibly affected by the disc gravity. We will also argue 
that the physical reason for this fact is the existence of three 
different time-scales involved in the problem: (i) the orbital 
period of the star around the central mass (short time-scale), 
(ii) the period of oscillations in eccentricity and inclination 
of the orbit (medium time-scale, these oscillations are due 
to the disc gravity), and (iii) the time for grinding the or- 
bital plane into the disc (long time-scale, due to successive 
interactions of the stars with the disc). Effects which can 
be ascribed to the medium time-scale present a new feature 
discussed in this paper within the context of galactic nu- 
clei surrounded by an accretion disc, although analogous ef- 
fects of oscillations or sudden changes in orbital parameters 
are well-known from other applications (cf. recent discus- 
sion on dynamics of planetary motion by Holman, Touma 
& Tremaine 1997; Lin & Ida 1997). 

Details of the model are described in the next section. 
Then, in Sec. ^ long-term evolution of stellar orbits is ex- 
amined. Finally, conclusions of our present paper are sum- 
marized in Sec. a. 



2 THE MODEL 

The discussion in the present paper is completely Newto- 
nian. Stellar orbits under consideration are bound to the 
central mass which determines their form but the orbits are 
not exact ellipses for two reasons which principally influence 
their long-term evolution, namely: 

(i) the gravitational field of the disc material acts on the 
star as a tiny perturbation; 

(ii) successive interactions with an accretion disc also affect 
trajectories in an impulsive manner, when the star intersects 
the disc plane. 

Each star is treated as a free test particle moving in the 
combined gravitational field of the centre and the disc; col- 
lisions with the disc material act as instantaneous periodic 
perturbations of the trajectory. All other effects are beyond 
the scope of the present paper (although they will have to 
be taken into account in a future self-consistent model). 

One can speculate that other effects, ignored in this pa- 
per, may have a comparable result on stellar orbits. In par- 
ticular, general relativistic dragging of inertial frames due 
to rotation of the central object and the disc material will 
break the spherical symmetry of the gravitational field and 
result in sudden excursions of the mean orbital parameters. 
On the other hand, gravitational radiation will result in a 
slow decay of the orbit in a manner analogous to the effect of 
direct collisions of the body of the star with the disc mate- 
rial. These effects also support our conclusion that star-disc 
collisions should not be considered as the only perturbation 
of stellar orbits when their long-term evolution is discussed. 
Nevertheless, our very restricted choice to the two effects 
listed above, (i) and (ii), has an additional reason. Indeed, 
these effects are linked to each other, because both of them 



are determined by the mass of the disc: increasing the mass 
of the disc affects the stellar trajectories more by its gravita- 
tional attraction, and, at the same time, star-disc collisions 
also become more important (being on average proportional 
to the surface density of the disc). We will demonstrate that 
a consistent model involving any one of the two effects must 
take the other effect into account too. But the main nov- 
elty of this paper is in even a stronger claim: it is the first 
of the effects mentioned above — the gravity of the disc 
— which influences long-term evolution of the stellar orbits 
dominantly, while collisions with the disc material represent 
an underlying mechanism causing a slow and continuous or- 
bital decay. In other words, changes in eccentricity and incli- 
nation are driven dominantly by the disc gravity and they 
can occur rather abruptly. From this perspective. Ranch's 
(1995) statistical model of the star-cluster evolution due to 
interactions with the central accretion disc requires also the 
inclusion of the influence of the disc gravity. 

It is also worth mentioning that, in analogy with Ranch 
(1995), we disregard, at this stage of the model, mutual 
interactions of the stars, considering the star cluster as a 
coUisionless system. A rigorous approach will require us to 
solve the Fokker-Planck equation in a manner analogous to 
Bahcall & Wolf (1976, 1977; see also Peebles 1972; Young 
1980; Shapiro & Teukolsky 1985, 1986; Zamir 1993; Quin- 
lan, Hernquist & Sigurdsson 1995; Sigurdsson, Hernquist & 
Quinlan 1995.) Obviously, this neglect represents a large 
simplification, especially in the nuclear region close to the 
central galactic object (Statler, Ostriker & Cohn 1987; Lee & 
Ostriker 1993). Nevertheless, demonstration of the effect of 
the disc gravity upon stellar orbits which we want to discuss 
hereafter does not call either for general theory of relativity 
nor mutual interactions among stars themselves to be taken 
into account. We expect, however, that the capture rate of 
stars by the central object can be only roughly estimated in 
this approach (the capture rate has been discussed in vari- 
ous approximation by, e.g., Frank & Rees 1976; Nolthenius 
& Katz 1982; Novikov, Pethick & Polnarev 1992; Hameury 
et al. 1994; Sigurdsson & Rees 1997). 

2.1 Gravitational field 

We describe the gravitational field as a superposition of the 
spherically symmetric field of the central object (potential 
Vc{r) — —GM^/r; r is radial distance from the centre) and 
an axially symmetric field of the disc (potential Vd{R,z); 
cylindrical coordinates {R,z), z = is the disc plane). The 
disc is geometrically thin and it is described by surface den- 
sity x{R). We assumed ^ — 10~^ for definiteness. It is 
worth mentioning that our >c{R) corresponds to vertically 
integrated density which is introduced in the standard the- 
ory of geometrically thin discs. With this correspondence 
one can compare our results with other works which employ 
Shakura-Syunaev (1973) and Novikov-Thorne (1973) discs. 
We do not consider a more complicated case of geometrically 
thick tori in this paper. 

An analytical expression for the disc potential can be 
found for some particular forms of the density distribution 
(Binney & Tremaine 1987, Evans & de Zeeuw 1992). Despite 
the fact that such models can approximate real astrophysi- 
cal discs only roughly, their main advantage is the analytical 
expressions for the gravitational field which they offer. Obvi- 
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ously, a careful check is needed in order to verify that none of 
the important qualitative features of the solution has been 
altered. We shall follow this line of reasoning by working 
mainly with highly simplified Kuzmin's class of discs. 

The surface density-potential pair for Kuzmin's model 
reads 



V4R,z) 



2n (A2 + i?2)3/2 ' 

GMd 



[R^ + {A + \z\f] 



21 1/2 



(1) 

(2) 



where A is a free parameter of the model. An easy exer- 
cise then yields components of acceleration due to the disc 
gravitational field. Kuzmin's discs are of infinite radial ex- 
tent but their mass is finite because the surface density de- 
creases with radius fast enough. Relevant formulae for discs 
of finite radial size and arbitrary surface-density profile are 
summarized in Appendix. 

In our numerical code for integration of stellar orbits we 
have used either the simple Kuzmin formulae given above or, 
in case of discs with a-priori unconstrained surface density 
distribution, pre-computed Vd and components of the gra- 
dient {dVd/dR,dVd/dz) for a given distribution of >f in a 
fixed grid of {R,z). Then we employed six-point interpola- 
tion formulae and evaluated the gravitational force acting 
on a star at any position. By performing several tests we 
have tuned parameters of the grid in order to optimize com- 
puter time necessary for integration of the stellar orbits and, 
simultaneously, to preserve a pre-determined accuracy. 



2.2 Interaction with the disc material 

The physics of collisions of the stars with the disc mate- 
rial can be represented by a prescription for the change of 
star's velocity. The prescription is based on a simplified hy- 
drodynamic scenario in which the star crossing the disc is 
treated as a body in hypersonic motion through fluid. The 
resulting change of velocity (which occurs always at z = 0, 
once or twice per each revolution) is obtained by integra- 
tion over a short period when star moves inside the disc 
material (as is done with all quantities in the thin-disc ap- 
proximation). Since the pioneering works of Hoyle & Lyt- 
tleton (1939), Chandrasekhar (1942), and Bondi & Hoyle 
(1944), the drag on a cosmic object has been considered for 
various astrophysical problems. It has been recognized that 
in the case of a supersonic flow which resembles our situa- 
tion (Ostriker 1983; Zurek, Siemiginowska & Colgate 1994, 
1996) the hydrodynamic drag consists of a component given 
directly by the fiow of material onto the stellar surface (or 
into a stellar-mass black hole; Petrich et al. 1989) and a long- 
range component given by the interaction of the flow with a 
conical or a bow shock surrounding the star. The relative im- 
portance of the two components is sensitive to the details of 
the flow as well as to complicated turbulent processes in the 
wake (e.g., Livio, Soker, Matsuda & Anzer 1991; Zurek et 
al. 1996). In what follows we will adopt a somewhat simpli- 
fied empirical model which, however, we argue still reflects, 
rather conservatively, physical effects in our consideration. 

We adopt a simplified formula for the mutual inter- 
action between star and disc as in Vokrouhlicky & Karas 
(1998). Conclusions drawn from the present paper can be 



thus compared with previous results in which the gravity of 
the disc was ignored (see also Syer et al. 1991; Artymowicz, 
Lin & Wampler 1993; Ranch 1995). The impulsive change 
of the star's velocity, 5v, when it crosses the disc will be 
described by 

5v ^T.{R,v)v,ri (3) 

with 'S{R,v) = ~Tv( >c {Rl/m*){vrei/v±), and C ~ 1 + 
(t;t/«)'' In A. Here, Vrei denotes the relative velocity of the 
star with respect to the disc matter, R is the radial coor- 
dinate in the disc, 7?* stellar radius, its mass and v-,, 
the escape velocity from surface of the star; v± is the nor- 
mal component of the star's velocity to the disc plane, and 
In A is the long-range interaction factor. Obviously, the lat- 
ter term is to be considered for transonic flows only, i.e. 
when the Mach number M ~ (R/H) > 1 (this condition is 
well-satisfied in standard thin accretion discs of Shakura & 
Sunyaev 1973; cf. also Zurek et al. 1994 who estimate the 
amount of the disc material swept out of the disc by a star's 
passages). Finally, H is the geometrical thickness of the ac- 
cretion disc at distance R. Ostriker (1983) gives a rough 
estimate of A « (H/_R*)(u/ii*)'^. Again, for standard thin 
discs one obtains A 3> 1. We note again that the main con- 
cern is not to underestimate the role of the hydrodynamic 
drag. We shall thus rather conservatively assume that the 
factor ^ is equal to 10"^. 

It is worth mentioning that Artymowicz (1994) proved 
formula (^) approximates also the interaction of the star 
with density and bending waves excited on the disc surface 
by the star itself (see also Hall, Clarke & Pringle 1996). This 
additional effect can be accommodated by an appropriate 
modification of the drag factor Taking into account an 
upper estimate on we effectively include this effect in our 
considerations too. 

To conclude this paragraph, we recall that the factor E 
in eq. (^ is proportional to the surface density and there- 
fore depends linearly on the total mass of the disc (given 
the density profile). When expressed in units of the central 
mass: E oc /i with a numerical factor depending uniquely on 
the characteristics of the moving object (neutron star, white 
dwarf, stripped star, etc.) and of the disc material. Though 
simplified, this model reflects the intuitive guess that more 
mass in the disc makes the effects of the hydrodynamic drag 
more profound. 



3 EVOLUTION OF STELLAR ORBITS 

Hereinafter, we examine stellar orbits. First, we discuss how 
the gravity of the disc influences individual trajectories (Sec. 
3.1). This will help us to illustrate the main differences with 
respect to previous works (especially. Ranch 1995). Next, 
we consider the combined effect of the disc gravity and the 
drag due to collision with the disc material (Sec. |j.2| ). In 
both cases, we will integrate orbits numerically and then 
we will describe orbital evolution in terms of osculating Ke- 
plerian elements. The elements relevant for our work are: 
semimajor axis a, eccentricity e, inclination / with respect 
to the disc plane (J > 90° corresponds to retrograde orbits 
while / < 90° corresponds to prograde orbits), and argu- 
ment of pericentre u as measured from the ascending node. 
The longitude of the node does not appear in the following 
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Figure 1. Contours C = const of the first integral ( |l2[ ) in tlie plane of non-singular variables (k, h). Radial distance from origin determines 
mean eccentricity e according to eqs. (13)— (14). Polar angle (measured from the horizontal fc-axis) is the argument of pericentre Q. A 
separatrix (thick curve) is drawn on the border between regions of different topology in panels (b)— (d); there is no separatrix in (a). 
Contour lines with polar angle acquiring values in the full range of < tD < 27r correspond to a zone of circulation, while some values 
of u) are not allowed in zones of libration. In this example, Kuzmin's disc with j4 = 75 (in units normalized by half of gravitational 
radius of the central mass) has been assumed. Graphs are shown for four values of Kozai's integral respresenting topologically different 
situations: (a) c = 0.9 (0 < e < 0.436); (b) c = 0.81 (0 < e < 0.586) — notice a figure of 8-shaped separatrix between the libration 
and circulation regions; (c) c = 0.7 (0 < e < 0.714) — notice the appearance of the second (inner) circulation region around the origin; 
(d) c = (0 < e < 1). In all these cases, trajectories have semimajor axis a = 200. Shaded area indicates highly eccentric orbits which 
must be trapped by the central object at some stage. More details are given in the text. 



discussion because of the axial symmetry of gravitational 
field. The (assumed) small value of the disc-mass parame- 
ter ^ and the value of C, guarantee that the time-scale of 
the orbital evolution is much longer than period of a sin- 
gle revolution. Therefore, we will average relevant quanti- 
ties over individual revolutions around the centre whenever 
it is appropriate. One can verify that /i controls the ratio of 
medium vs. short time-scale while C, affects the ratio of long 
vs. medium time-scale. We accept physically substantiated 
values for ^ (< 10"'^) and (1-10^) which guarantee that 
the three time-scales are well-separated from each other. 

3.1 Effects of the disc gravity 

The motion of cosmic objects in gravitational fields of discs 
or rings of matter has been considered in the context of solar 
system studies (e.g., Ward 1981; Heisler & Tremaine 1986; 
Lemaitre & Dubru 1991; McKinnon & Leith 1995), recently 



discovered planetary systems (Holman, Touma & Tremaine 
1997), and in galactic dynamics (Huang & Carlberg 1997). 

As mentioned above, we are interested in the long-term 
variations of orbits around a massive centre and a much less 
massive disc. The averaging technique is a very useful ap- 
proach for understanding the qualitative behaviour of the 
system (Arnold 1989). It can be formalized in terms of a se- 
ries of successive canonical transformations in which rapidly 
changing variables (e.g. mean anomaly along the osculating 
ellipse) are eliminated (Brouwer & Clemence 1961). Since 
the reader may not be closely familiar with this technique 
we briefly remind several basic concepts. 

Suppose we consider motion of a particle in the Ke- 
plerian fleld, determined by an integrable Hamiltonian Hq 
and a weak perturbation. Perturbation is described by a 
potential eV . The parameter e indicates a smallness of the 
perturbation (in our case, the disc/hole mass ratio /j plays 
the role of e). Hamiltonian of the central field in terms of 
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Figure 2. The same as in Figure y but for orbits with a smaller 
(b) c = 0.84; (c) c = 0.7; (d) c = 0. 

the Keplerian elements reads Ho — —GMc/{2a), depending 
uniquely on the semimajor axis, while the perturbing po- 
tential V depends, in general, on all Keplerian parameters. 
For this reason, a general solution cannot be found but the 
the problem is simplified when some symmetry appears. For 
instance, the axial symmetry induces independence on the 
longitude of node. Since the averaging technique is based on 
canonical transformation theory one should rather use a set 
of canonical elements of the two-body problem. Delaunay 
parameters are the most common choice (see, e.g., Brouwer 
& Clemence 1961). 

In order to grasp the long-term evolution of the system 
one needs to get rid of the "fast" (i.e. rapidly changing) vari- 
ables. In our problem there is the only one fast variable: the 
mean anomaly. The idea of the averaging method is then 
formalized by seeking a set of new Delaunay variables using 
a canonical transformation under the constrain of indepen- 
dence of the transformed Hamiltonian on this fast variable. 
It is possible, indeed, by formal development in the small pa- 
rameter e. The original perturbing potential V reads, after 
the transformation, 

eV ^ eV = eVi + eV2 + e^Va + ■ ■ • . (4) 

Various iterative and computer-algebra adapted methods for 
recursive generation of the potentials Vi, V2, etc. were de- 



(b) 




k 



semimajor axis, a = 75. Values of Kozai's integral are: (a) c = 0.9; 



veloped (e.g. Hori 1966; Deprit 1969). We note that the first 
term in this series, notably Vi, is just the average of the 
original potential V over the mean anomaly. Obviously, a 
rigorous proof of convergence of the series (W), and thus a 
success of the whole procedure, remains a very difficult task. 
Despite of this fact the averaging technique is often very use- 
ful. Typically, results based even on highly truncated part of 
this series are valid on a limited time interval though they 
fail to predict the system's evolution on an infinite time- 
scale. In many applications such a weakened demand is suf- 
ficient (see, for instance, Milani & Knezevic 1991 for the 
application to a long-term evolution of asteroidal orbits). 

Regarding the independence of the new Hamiltonian 
on the fast variable (i.e. mean anomaly) we immediately 
conclude that the mean perturbing potential V is the first 
integral of motion: 

eV ^eVi+e^V2+e''V3 + ...^eCie) . (5) 

Factorizing out the small parameter e one can see that the 
first order perturbing function, Vi , is nearly constant — pro- 
vided that the higher order contribution is neglected. Notice 
that the right-hand side constant C is a function of the small 
parameter e. However, in the most truncated level of aver- 
aging, accounting for the first order term of the right hand 
side of (^) only, we have 
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Figure 3. Mean eccentricity e vs. semimajor axis d. The curly 
curv e co rresponds to the trajectory from the first example in 
sec. 3.2. Both axes are logarithmic. Kuzmin's disc model with 
scaling parameter A = 75 and mass parameter /i = 10~^ has 
been used. Initial semimajor axis d = 200, eccentricity e = 0.6, 
inclination / = 34.4° and argument of pericentre uj = 90? For sake 
of comparison, three hypothetical orbits, evolution of which disre- 
gards gravitational influence of the disc, are also plotted (mono- 
tonic curves); one of these curves (labeled by "s") corresponds to 
the same initial conditions as the complete solution. 



Vi ^ C(0) + 0{e) 



(6) 



and the averaged perturbation potential Vi is (approxi- 
mately) constant. In the following we shall drop out the 
argument (which is zero) and write C for simplicity. In the 
case of an axially symmetric system the integral (^) is suffi- 
cient for global integrability. Now the problem is reduced to 
a single degree of freedom and the whole phase space can be 
plotted in a simple two-dimensional graph. Qualitative fea- 
tures of the solution can be illustrated by isocurves of the 
C-integral (^. 

In the remaining part of this section we shall apply the 
previous brief review of the averaging to our problem. We 
shall restrict ourselves to the first-order procedure in which 
the perturbation potential Vd is substituted by its average 
Vd over the mean anomaly: 



1 



dv 



Vi{R,z) 



(7) 



where r, R, and z are functions of true anomaly v, rj = 
^\ — e? (e denotes mean eccentricity of the orbit). One can 
write, in terms of mean inclination / and mean argument of 
pericentre w, 



1 -f e cos V 
z = r sin I sin {p -\-v) , 

R = r\/ 1 — sin^ Jsin^ {uj^v). 



(8) 

(9) 
(10) 



At this level of approximation, the mean semimajor axis 
a stays constant and differs from the osculating semimajor 



axis o by short-period terms only. Owing to axial symmetry 
of Vd there exists an additional first integral of motion which 
relates I to the corresponding value of e in the course of their 
evolution, 



yj 1 — cos / : 



; const . 



(11) 



Eq. ( [li[ ) is often called Kozai's integral (Kozai 1962). Here 
again, the overbar distinguishes the mean elements from the 
corresponding osculating elements. The problem is reduced 
to the evolution of e and uj which are constrained by 



Vd (e,tj; c,a) = C . 



(12) 



We will introduce a pair of non-singular canonical variables 
{k,h): 



k = 



h 



y^2(i- yr^^' 



(13) 
(14) 



(fc, h) are often called Poincare variables, since they have 
been introduced to orbital dynamics by Poincare (1892). 
Levels of C = const in the (k, /i)-plane offer a convenient 
representation of the long-term evolution of mean orbital el- 
ements. At small values of eccentricity e, the radial distance 
from the origin in the (fc, /i)-plane is equal to the mean ec- 
centricity itself, while the polar angle has the meaning of the 
argument of pericentre ui. The final task is to evaluate the 
averaged potential, Vd (e, tli;c, a). Since we are interested in 
averaging over orbits which intersect the disc, no simple an- 
alytical techniques based on expansion in orbital elements 
(common in celestial mechanics) can be applied to a gen- 
eral distribution of surface density >c{R) (see the Appendix). 
Even in the case of simple Kuzmin's discs (2) the averaging 
cannot be performed in analytical functions and must be 
obtained numerically. 

Figure ^ illustrates contours C = const of the first inte- 
gral ( |l2[ ) in the case of a Kuzmin disc with scaling parameter 
A = 75. The size of orbit is characterized by a value of mean 
semimajor axis a = 200. Geometrized units of length have 
been used; in the following, we will use half of the gravita- 
tional radius of the central object as a natural unit of length 
(this choice is motivated by interpreting the central object 
as a black hole; otherwise, our discussion is Newtonian). 

Kozai's integral c spans the whole interval (—1, 1), but 
we can restrict ourselves to positive values because c occurs 
in Vd only squared. Large values of c (Fig. |^a, c = 0.9) corre- 
spond to quasi-circular orbits with low inclination to the disc 
plane and small oscillations in mean eccentricity. The argu- 
ment of pericentre circulates in the whole interval (0,27r). 
When c is decreased (Fig. ^o, c = 0.7), larger inclinations 
occur and two regions of libration with associated stable 
points at a) = ±7r/2 develop and bifurcate into a figure of 
8-shaped region. Orbits outside this libration region still cir- 
culate with < a) < 27r (0 is to be identified with 27r). Notice 
that circulating trajectories which are close to the separa- 
trix of the two regions exhibit large oscillations in the mean 
eccentricity. At still smaller values of c (Fig. |l|c, c = 0.2), 
the circulation region bifurcates in its inner part near origin. 
Stable points in the libration regions are expelled farther to 
higher eccentricities. Setting c = (Fig. |l|d) the inclination 
is constrained to I — n/2 (polar orbits). There is a maxi- 
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mum eccentricity above which the orbits are trapped by the 
central object (this situation corresponds to a pericentre dis- 
tance equal to 2 in our example). We observe that a large 
portion of the plot (indicated by shading) corresponds to or- 
bits which emerge from or fall into the centre, while only a 
very minor part — namely the inner circulation region with 
a stable point e = — contains orbits which survive the 
long-term evolution. In the following, this property turns 
out to be essential for the fate of orbits counter-rotating 
with respect to the disc. We will argue that most initially 
retrograde orbits cannot last long enough to be tilted over 
the polar orbit and inclined into the disc plane because, in 
course of this process, their eccentricity is being pumped up 
to such a high value that they get captured by the central 
object. 

Figure |^ shows the same graphs as in Figure |l] but for 
different ratio of the scaling parameters of density distribu- 
tion {A = 75, as before) and the orbit (a = 75, changed). All 
features in the (fe, /i)-plane described before stay unchanged 
except the fact that libration regions and the inner zone of 
circulation cover a larger portion of the graphs. At small val- 
ues of Kozai's integral we again observe that eccentricity is 
forced to large oscillations which eventually lead to capture 
(shaded area). 

We have verified that qualitatively similar results hold 
for different, astrophysically relevant disc density-potential 
pairs. In particular, using formulae from the Appendix we 
have numerically computed potential of discs with -R~^, 
and surface density laws and then averaged over corre- 
sponding quasi-elliptic orbits. 

3.2 Influence of star-disc collisions 

We now consider inclusion in our model of the influence of 
the hydrodynamic drag affecting the orbit when the star 
crosses the accretion disc. We will use the approximations 
described in Sec. 



2.2 



and, for definiteness, we will assume 
that the central cluster of stars is formed by white dwarfs 
or stripped cores of ordinary stars and use relevant param- 
eters for estimating numerical factor in the relation E oc /i 
(eq. above). Before embarking on a description of our re- 
sults, we recall that Ranch (1995) has considered long-term 
evolution of stellar orbits in which he took into account the 
impulsive drag acting twice per revolution only (no gravity 
of the disc; see also previous studies by Syer et al. 1991, and 
Vokrouhlicky & Karas 1993, 1998). 

All stellar orbits interacting with the disc, indepen- 
dently of their initial conditions, exhibit long-term decay of 
the semimajor axis, and monotonous decrease of eccentric- 
ity and inclination (grinding to the disc plane). For initially 
low-inclination orbits, the characteristic time-scale of circu- 
larization is comparable to the grinding time after which 
the orbital plane becomes inclined into the disc. Orbits with 
large initial inclination, in particular all initially retrograde 
orbits, are circularized before they incline to the disc. Ranch 
(1995) argues that, for orbits with an initially moderate ec- 
centricity, a particular combination of the mean elements, 
a(l — e^) cos*(//2), is quasi-conserved in the course of evolu- 
tion (see also Vokrouhlicky & Karas 1998). This property is 
insensitive to a particular model of the star-disc interaction, 
surface density profile, and the total mass of the disc. For 
further use we recall the way that the hydrodynamical drag 



affects the mean semimajor axis a and the Kozai parameter 
c (two integrals of the long-term evolution when star-disc in- 
teractions are neglected): due to the drag effect, a undergoes 
permanent decay while c typically increases from its initial 
value to the final value of Cf = 1, corresponding to a circular 
orbit in the disc plane. Initially retrograde orbits behave in 
a somewhat different way: their circularization time-scale is 
significantly shorter than the grinding time. For these or- 
bits, the value of Kozai's parameter first slightly decreases 
before exhibiting monotonous increase to Cf. 

The main result of the present paper is that most of the 
above-mentioned properties cease to he true in our general- 
ized model which, besides star-disc collisions, includes also 
the gravitational influence of the disc matter. The key point 
is the fact that the characteristic time-scale a/a ~ c/c of the 
hydrodynamical-drag effects is much longer than the charac- 
teristic time-scale for circulation and libration of pericentre. 
Circulation and libration along C = const lines are due to 
the disc gravity, and the corresponding period is also the 
time-scale for gravity-driven oscillations in eccentricity and 
inclination. As a consequence, the principal features of or- 
bital variations explained in the previous section remain un- 
changed apart from a very slow adiabatic evolution of quasi- 
integrals (a, c). However slow this evolution is, stellar orbits 
are strongly and abruptly affected at some stages: when they 
cross the separatrix (due to some perturbation) or when the 
libration region bifurcates. In the following paragraphs we 
will illustrate these facts by showing typical orbits. Kuzmin's 
disc with A = 75 has been assumed in examples described 
below. The mass of the disc has been set to fj, = 10~^ for def- 
initeness (units of the central mass). It should be mentioned 
that the results depend only very weakly on the particu- 
lar value of fi, provided it is sufficiently small. The reason 
emerges from the fact that the disc-mass parameter can be 
factorized out of all expressions containing the averaged po- 
tential Vd (which determines all important dynamical fea- 
tures). On the other hand, and in contrast to the results of 
Ranch (1995), one expects dependence on the surface den- 
sity profile k{R). 

In the first example of this section, we consider an orbit 
with the initial parameters: d — 200 (in units of one half of 
the gravitational radius), e = 0.6, I = 34.4° (c « 0.66), and 
uj — 90? Figure |^ shows mean eccentricity e as a function 
of the semimajor axis d during the course of orbital evo- 
lution. For sake of comparison we have also plotted three 
curves of orbital evolution if the gravitational influence of 
the disc is ignored. These curves can be compared directly 
with corresponding results of Vokrouhlicky & Karas (1993) 
and Ranch (1995). The curve labeled by "s" corresponds 
to the same initial conditions as the orbit of the complete 
model with gravitational effects taken into account (notice 
the difference in predicted radius of the terminal orbit); the 
other two curves correspond to different initial eccentric- 
ities. As expected from the previous discussion, the most 
evident features are (i) large oscillations in eccentricity, and 
(ii) permanent decay of the semimajor axis. The radius of 
the resulting circularized orbit in the disc plane, flf = 120.2, 
differs slightly from the estimation based on Ranch's (1995) 
formula, af' ~ ai(l — ef)cos'*/i/2 — 106.5. A closer look at 
orbital evolution explains a remarkable bump in eccentric- 
ity near d ~ 141. In the initial state, the pericentre of the 
orbit librates in one of the two lobes of the 8-shaped zone. 
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Figure 4. Long-term evolution of the orbit from Figure js] projected onto the {k, /i)-plane. Initial librations around the fixed point at 
u> = 90° are shown in panel (a). As the Kozai's parameter c increases, the libration zone shrinks and the orbit approaches the separatrix. 
Eventually, the orbit crosses the separatrix and starts circulating around the origin; see continuation of orbital evolution in panel (b). 
Still further continuation in (c) shows that the orbit terminates in the inner zone of circulation, with zero eccentricity. 



around Co = 90° as shown in Figure As the Kozai param- 
eter c slowly increases due to the hydrodynamic drag the 
libration zone shrinks and the orbit approaches the separa- 
trix. At a particular instant, the orbit is expelled from the 
libration zone, crosses the separatrix and enters the zone of 
circulation (see Figure ^). This transition is accompanied 
by large oscillations of the mean eccentricity and temporal 
slowing down of the eccentricity decay. 

Our second example of the orbital evolution (Figure ^ 
appears even more peculiar (from the viewpoint of previous 
results neglecting gravity of the disc) . Initial orbital param- 
eters are as follows: a = 200 , e = 0.638, I = 67° (c « 0.3), 
and w = 0. The first part of the orbital evolution ends at 
a ~ 71.5 by a significant increase of mean eccentricity. The 
oscillations of eccentricity eventually settle down and the 
evolution terminates as a circular orbit in the disc plane, 
ttf ~ 49.2. The latter radius is to be compared with Ranch's 
approximative result: af" « 54.3. Obviously, this estimate 
fails to predict the final radius correctly but the difference is 
still within uncertainty of Ranch's quasi-integral (Vokrouh- 



licky & Karas 1998). Again, a closer look at the evolution 
of pericentre in the (fc, /i)-plane sheds light on properties of 
this model. Initially the orbit is confined to the inner zone 
of circulation around the origin (Figure |^ . An adiabatic in- 
crease of Kozai's parameter results in collapse of this zone 
and the orbit is expelled towards the separatrix. At the in- 
stant of crossing the separatrix, the orbit may either enter 
the libration zone surviving for larger values of c, or it ter- 
minates in the large circulation zone. It can be argued, that 
orbits close to the separatrix spend most of their time near 
the hyperbolic points which leads preferentially to the cap- 
ture in the circulation zone. This indeed happened in our 
example, as indicated by a thick line in Figure ^ The form 
of the 8-shaped libration zone results in significant increase 
of eccentricity. Figure |^ shows a projection of the orbit onto 
the plane of mean eccentricity e and inclination I. One can 
notice that individual oscillations are confined to the under- 
lying grid of constant Kozai's integral (pT|). Slow (adiabatic) 
diffusion across the lines c = const reflects a long-term fea- 
ture of the orbital evolution. Transition from the inner to the 
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Figure 6. Orbital evolution from Figure ^projected onto the (fc, /i)-plane. Initial circulation in the inner zone of circulation (thin curve) 
is followed by transition to the external zone of circulation (thick curve) at a critical value c fs: 0.78 of the Kozai parameter. At this 
moment eccentricity is increased back to a high value. Notice also a reversal of the sense of circulation. 



outer circulation zones is accompanied by a large increase 
of eccentricity, but only a moderate decrease of inclination. 

So far we have demonstrated the intricate role of the 
disc gravity in evolution of stellar orbits with initially pro- 
grade inclination (Ji < 90°). Even though details of the or- 
bital evolution are different when compared to models ne- 
glecting disc gravity, overall features remain approximately 
unchanged. Most importantly, radii of circularized orbits in 
the disc plane are comparable. However, as mentioned be- 
fore, one has to pay particular attention to orbits with ini- 
tially large, retrograde inclination where the results become 
qualitatively and quantitatively different. Two subsequent 
examples deal with this class of orbits. 

Figure ^ corresponds to an initially retrograde orbit 
with parameters: a = 200, e = 0.1, / = 120° (c ^ -0.5), 
and ijj = 90? Similar to the previous example, this orbit is 
originally locked in the inner zone of circulation. This is nec- 
essary (but insufficient as we shall see later) if the orbit is to 
be tilted over the polar orbit to a prograde one and survive 
further evolution. Otherwise, it gets captured rather soon 
by the central mass. Just before the separatrix of the inner 



circulation zone shrinks to origin, the orbit is released to the 
outer zone of circulation. The existence of the libration lobes 
then leads to significant increase of the mean eccentricity, up 
to 0.9. The (fc, /i)-plane representation of the orbit evolution 
is given in Figure |^. The circular equatorial orbit with ra- 
dius Of ~ 7.8 is a final state of this evolution. Its radius is 
to be compared with af" ~ 12.4 (no gravity of the disc). The 
difference between terminal radii predicted by the two mod- 
els amounts to 30 %. Starting with various initial conditions 
we found that terminal radius and corresponding time are 
comparable, typically within factor of 2. 

A truly fundamental difference between the complete 
model, with the disc gravity taken into account, and the 
simplified model, disregarding effects of the disc gravity, is 
observed when the initial eccentricity of the previous orbit 
is slightly increased to 0.3. The corresponding orbital evolu- 
tion is shown in Figure |l^. The orbit is again initially locked 
in the inner zone of circulation. But now, at the transition 
to the outer circulation zone, eccentricity increases over a 
critical level and the orbit is captured by the central object 
(pericentre is less than 2). On the contrary, a hypothetical 
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Figure 5. Eccentricity e vs. semimajor axis 5 (the second exam- 
ple in sec. |3.2| ). The same parameters of the disc as in Figure ^ 
different initial parameters of the orbit: a = 200, e = 0.67, 1 = 67° 
and (Ii = 0. Again, the gravitation of the disc induces oscillations 
and an abrupt increase of eccentricity (curly curve) which can- 
not be seen when the disc gravity is neglected (two monotonic 
curves). 

orbit with no gravity of the disc grinds to the disc plane 
at final radius Of « 11.4. Concerning retrograde orbits, we 
can conclude that even those ones that have started within 
the inner zone of circulation are not safe from capture (ob- 
viously, orbits with larger initial eccentricity than 0.3 in our 
example are also captured) . We have also verified that retro- 
grade orbits which are initially locked in the libration lobes 
or in the outer circulation zone do not survive tilting over 
the polar orbit, being soon captured by the centre due to 
large oscillations in eccentricity. All these orbits are grinded 
safely to the equatorial plane in the framework of the sim- 
plified model of Rauch (1995), when the disc gravity is ne- 
glected. We would like to stress again that this difference is 
not based on the particular value of the disc mass we have 
chosen in our examples (/i — 10"''); rather it is present for 
an arbitrary nonzero mass of the disc (naturally, evolution 
takes place on a longer time-scale for smaller values of fj). 
The same results hold also for less massive discs. 

A theoretical substantiation for the reported /.i- 
independence of our results is based on existence of the 
quasi-integral (^ mentioned above. Notice that the small 
parameter e (i.e. /i in our case) can be factorized from this 
formula. Only when the higher-order terms of the exact inte- 
gral (^ are taken into account the results (phase portraits, 
in particular) depend on fj,. As soon as ji is sufficiently 
small, which means smaller than about 10~^ in practice, 
the averaging technique offers a simple explanation for the 
^-insensitivity of the results. 

We have verified our principal conclusions for discs with 
different surface density profiles. In particular, we considered 
stellar orbits around discs with density profiles proportional 
to , R~^, and R~^. Obviously, the values of parame- 
ters when transitions between zones in the (fc, /i)-plane oc- 
curred are quantitatively different in individual cases, but 
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Figure 7. Mean inclination I vs. mean eccentricity e of the orbit 
from Figure ^ Background grid of constant values of the Kozai 
parameter, c = const, is also plotted. Orbital evolution consists 
of fast oscillations along the lines of c = const and slow diffusion 
towards the maximum value of Cf = 1 (e = J = 0). Kozai's 
quasi-integral is not isolating because the transition between the 
inner and the outer (shaded) circulation zones crosses one of the 
c = const curves. Rauch's solution is also shown by the monotonic 
curve (labeled by "R"). 
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Figure 8. Mean eccentricity e vs. semimajor axis a (the third 
example in the text; an initially retrograde orbit). The same pa- 
rameters of the disc as in Figure jsj starting parameters of the 
orbit are as follows: a = 200, e = 0.1, / = 120° and O = 90? 
Complete model orbit is represented by a curly curve with oscil- 
lations in eccentricity, while the two orbits of a simplified model 
(no gravity of the disc) exhibit monotonic decrease of eccentricity 
during the whole evolution. 
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what holds unchanged are the quahtative results. Most im- 
portantly, a significant fraction of the initially retrograde 
orbits is captured by the central mass. 

3.3 Notes on statistical properties of the cluster 

In this section, we briefly comment on possible influence of 
our results on statistical properties of a cluster of stars near 
a galactic nucleus. 

The important finding, in this respect, concerns the sig- 
nificant portion of orbits which increase eccentricity and 
than become captured by the centre, many of them being 
initially retrograde. These orbits arc missing in the final con- 
figuration of the system. The form of this final configuration 
appears to be sensitive to the detailed nature of the situar 
tion under consideration. Ranch (1995) specified the initial 
distribution of stars as a power-law in binding energy. Then 
the system evolved under interactions with a disc for infi- 
nite time, stellar orbits were inclined into the disc plane or 
captured by the centre, and the index of the power-law has 
been evolved accordingly. No new stars have been inserted 
into the system. Since retrograde orbits end up with smaller 
radii in the disc plane than initially prograde ones [due to 
the term cos"* (7/2) in Ranch's estimate], we expect in our 
model deficiency of orbits very close to the central object 
when compared to the results of Rauch. As a consequence a 
softer power law of the final distribution is to be expected 
in comparison to the Rauch's work. 

Alternatively, one can look for equilibirum in which the 
number density of stars remains constant at all radii. Stars 
which get captured by the central mass must be substituted 
by new ones which come from infinite radius. Now it is im- 
portant to take into account the fact that time interval for 
grinding the orbit into the disc depends on initial inclination 
(it is longer for initially retrograde orbits than for prograde). 
Discussion in this section is to be verified by detailed numer- 
ical simulations (work in progress). 



4 CONCLUSIONS 

It has been recognized by previous works that statistical 
properties of central galactic clusters are influenced by an 
accretion disc surrounding the nucleus because of twice-per- 
revolution interaction which affects stellar motion. The main 
results of this paper can be summarized as follows: 

(i) we demonstrated that any consistent model of the star- 
disc interaction has to take the influence of the disc gravity 
into account, in addition to the effects of direct collisions 
with gaseous material; 

(ii) as a result of the disc gravity, individual stellar orbits 
exhibit evolution which is different if compared to the situ- 
ation when collisions are considered but gravity neglected. 
Most importantly, we found that a significant fraction of 
initially retrograde (i.e. counter-rotating with respect to the 
disc material) orbits are captured by the central object. This 
is duo to large oscillations in eccentricity which affect polar 
orbits. 

We wish to note that our two findings mentioned above 
are to some extent different in their nature. The former one 
— (i) — is essentially a statement of consistency claiming 



that any reasonable model which involves the effects of di- 
rect star-disc physical interaction has to take disc gravity 
also into account. We argued that this claim is valid for all 
astrophysically reasonable objects expected in central clus- 
ters of galaxies: neutron stars, white dwarfs and stripped 
stars. The logic behind this result is due to the fact that 
both effects are controlled by the total mass of the disc. 
The latter finding — (ii) — then states how the model sup- 
plemented by effects of the disc gravity differs from previ- 
ous simpler models. Gravity of the disc induces dynamical 
structures, libration and circulation zones of the argument 
of pericentre. Adiabatic change of quasi-integral quantities 
and related transitions of trajectories between the two zones 
are the essence of our results. 

It is worth recalling that the above-mentioned results 
did show sensitivity on a particular model of the disc, espe- 
cially on the radial gradient of the surface density. Indeed, 
while Rauch (1995), considering only star-disc collisions, re- 
ported his results to be insensitive to a particular profile of 
the surface density or even to the model of the star-disc in- 
teraction, we observed that the fraction of retrograde orbits 
captured by the central mass in the course of their evolution 
depends on details of both star-disc collisions and effects of 
the disc gravity. On the other hand, our results show only 
a weak sensitivity on the total mass of the accretion disc. 
This feature can be easily understood by realizing that the 
disc mass parameter jj, factorizes out (in the first order of 
approximation) from the averaged potential Vd. As a con- 
sequence, the value of /i = 10~^ taken in our examples in 
Sec. 3 is not essential, and similar results hold also for less 
massive discs. 

Note: we have prepared a Java animation which illus- 
trates long-term evolution of stellar orbits in the two zones 
of the (&, /i)-plane (libration and oscillation in eccentricity); 
cf. "http://astro.troja.mff.cuni.cz/karas/papers/discapplet. 
html". 
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Figure 9. Long-term evolution of tlie orbit from Figure g projected onto the (fc, /i)-plane. Evolution in the inner zone of circulation 
is shown in panel (a), until the trajectory escapes to the outer zone of circulation. Form of the 8-shaped separatrix is indicated at 
the moment of transition to high eccentricity. Subsequent evolution continues in panel (b) with a steady decrease of eccentricity. The 
trajectory eventually approaches the origin of the [k, h) plane. 
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APPENDIX A: GRAVITATIONAL FIELD OF A 
THIN DISC 

In this section we describe a general method for evaluating 
the gravitational potential, and its gradient, of an axisym- 
metric disc. We introduce cylindrical coordinates {R, z, </>) 
with origin in the centre of the disc and the plane z = 
coinciding with the disc plane. The gravitational potential 
of the disc (outer radius fed) evaluated at arbitrary position 
{R,z) is given by (see, e.g., Binney & Tremaine 1987) 

V,{R,z)^-4Gj^'^^^^K[k{R')]dR', (Al) 
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Figure 10. Ecce ntricity e vs. semimajor axis a (the fourth ex- 
ample in sec. |3.2| ; initially retrograde orbit). The same disc and 
orbit parameters as in Figure ^ except for the initial mean eccen- 
tricity e = 0.3. The orbit which follows from our complete model 
has been captured by the central mass due to a large increase 
of eccentricity at the moment of transition from the inner to the 
outer zone of circulation, while the simplified model orbit grinded 
down safely to the equatorial plane. 



where 
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and 
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Here, K{k) is the elliptic integral of the first kind. (The 
dependence of functions B{R') and k{R') on coordinates R 
and z is not indicated explicitly.) 

No apparent reduction of integral ( |A1[ ) is possible until 
the radial distribution of the density >c{R') is specified. An 
important example which we will need below is a uniform 
disc, x{R') = Xo- The potential (Al) expressed in terms of 
elliptic integrals K, E and D reads, for R < bd, (Lass & 
Blitzer 1983) 
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Formula (A4) holds also in the region R > bd provided the 
first term in brackets of the right-hand side is suppressed. 
The expression for the potential inside the disc {z = Q, R < 
bd) can be simplified by applying the Gauss transformation 
of elliptic functions (Byrd & Friedman 1971). One finds that 



a more compact formula than the one given by Lass & 
Blitzer (1983). 

The components of grav itational force are given by the 
gradient of the potential ( Al). Direct algebraic manipulation 
results in 



dVd 
dR 



dVd 

dz 



2G ^<mKiE\kiR!M?^l^^±^ 
rI B[R') I^LM^jJ AHR') 



(A7) 
(A8) 

(A9) 



-K [k [r!)] jdi?', 

><{R')R'E[k{R')] 
-^^^X B{R')AHR') ''''' 



where we denoted 

{R')^z^ + {R^Rf . 
For a uniform disc one obtains 



dR 



dz 



2Gxq 



RBibdA'^'^''^'^'^^''^'^^^ 
{z^ + bl + R^)K[k{bd)]Y 



2GXo I ± TT 

bd~R 



+ 



bd + R 



n [q' {bd) ; k {bd)] 



(AlO) 



(All) 



The integrands in eqs. (Al) and (A7)-(A8) diverge in the 
disc plane, z ^ Q, although the result of integration must 
be finite because the potential is continuous across the disc. 
Taking into account relation 



K{k) ocln (l - k^) 



(A12) 



for fc ~ 1, one concludes the divergence in the poten- 
tial is proportional \nz. To get rid of numerical errors one 
conveniently splits the integrals into two parts by setting 
x{R') = [x{R') - k{R)] + k{R). Then the potential is 

Vd{R,z) = -4Gj^'^l^<2^^^^^WLK[k{R')]dR' 



+Vu [R, z; >c{R)] , 



(A13) 



where the second term corresponds to the potential of the 
disc with constant density >c{R) given by eq. (A4). Now 
the integrand in (A13) is well behaved. A sharp decrease 
of the integrand near R' ^ R, z m can be treated by 
appropriate numerical methods. The same approach can be 
applied successfully to evaluate the components of force in 
eqs. (A7)-(A8). 

Previous formulae are easily generalized to the case 
when the inner edge of t he d isc is at radius fld 7^ 0. The 
lower integration limit in ( |ai| ) and (A7)-(A8) is changed to 
ad, and in the case of formulae (A4) and (AlO)-(All) for 
the uniform density disc one employs a superposition of the 
field of a fictitious uniform disc with radius Od and formal 
density —Xo- 
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